#include "dat3d.h"

double angulito(double alfa0, double delta0, double alfa1,double delta1)
{
        double alfaMR;
        alfaMR=sin(delta0*M_PI/180)*sin(delta1/180*M_PI)
                +cos(delta0*M_PI/180)*cos(delta1/180*M_PI)*
                cos(M_PI*(alfa1-alfa0)/180);///OJO con las unidades!! de alfa
        if(alfaMR>1)alfaMR=1;
        else if(alfaMR<-1) alfaMR=-1;
        return (double)180*(acos(alfaMR))/M_PI;

};
double angul(double alfa0, double delta0, double alfa1,double delta1)
{
        double alfaMR;
        alfaMR=sin(delta0*M_PI/180)*sin(delta1/180*M_PI)
                +cos(delta0*M_PI/180)*cos(delta1/180*M_PI)*
                cos(M_PI*(alfa1-alfa0)/180);///OJO con las unidades!! de alfa
        if(alfaMR>1)alfaMR=1;
        else if(alfaMR<-1) alfaMR=-1;
        return (double)(acos(alfaMR));

};
double angul_z_f(double theta, double phi, double theta1,double phi1)
{
        double alfaMR;
        alfaMR=cos(theta)*cos(theta1)
                +sin(theta)*sin(theta1)*
                cos(phi-phi1);///OJO con las unidades!! de TODO
        if(alfaMR>1)alfaMR=1;
        else if(alfaMR<-1) alfaMR=-1;
        return (double)(acos(alfaMR));

};

double angul_point(pointing A, pointing B)
{
        double alfaMR;
        alfaMR=cos(A.theta)*cos(B.theta)
                +sin(A.theta)*sin(B.theta)*
                cos(A.phi-B.phi);///OJO con las unidades!! de TODO
        if(alfaMR>1)alfaMR=1;
        else if(alfaMR<-1) alfaMR=-1;
        return (double)(acos(alfaMR));

};
double align_point(pointing A, double dis,double ap,pointing B)
{
        double x,y,res;
	double del;
	del=A.phi-B.phi;
	x=((cos(B.theta)-cos(A.theta)*cos(dis))/(sin(A.theta)*sin(dis)));	
	y=(sin(B.theta)*sin(del)/sin(dis));
	res=atan2(y,x);
/*	if(del>0)res=acos(x);
	else res=M_PI-acos(x);*/
	if(res<0)res=M_PI+res;
	res=fabs(res-ap);
	if(res>M_PI/2) res=M_PI-res;
        return (double)res;
};

long double fchiqui(char separador,FILE * arch)
{
        char ch[600];
        int j;
	for(j=0;!feof(arch) && ch[0]==' ';ch[0]=fgetc(arch));
        for(;!feof(arch) && ch[j-1]!=separador;ch[j++]=fgetc(arch));
        ch[j-1]='\0';
//	printf("%i:%s\n",j,ch);
        return strtold(ch,(char **) NULL);
};
void grafica(char * NOM,int paso,int tot)
{
        div_t ppp;
        char g[4];
	g[0]='/';
	g[1]='-';
	g[2]='\\';
	g[3]='|';

        ppp=div(paso,4);
        printf("\r%s%c..........................%c %3i%%",NOM
        ,g[ppp.rem],g[ppp.rem],(int)((paso*100/(tot))));
};
void graf2(char * NOM,int paso,int tot,int NN=26)
{
        div_t ppp;
	int i;	    
        char g[4];
	g[0]='/';
	g[1]='-';
	g[2]='\\';
	g[3]='|';

        ppp=div(paso,4);
		printf("\r%s%c",NOM,g[ppp.rem]);
		for(i=0;i<NN;i++){
		if(i<paso*NN/tot) putchar('.');
		else if(i==(int)ceil((double)paso*NN/tot)) putchar('>');
		else putchar(' ');
	}
	printf("%3i%%",(int)paso*100/tot);
};
long double fchiqui2(char separador,FILE * arch)
{
        char ch[600];
        int j;

//	for(;!feof(arch) && ch[0]==' ';ch[0]=fgetc(arch));
        for(j=0;!feof(arch) && ch[j-1]!='.';ch[j++]=fgetc(arch));
        for(;!feof(arch) && ch[j-1]!=separador;ch[j++]=fgetc(arch));
        ch[j-1]='\0';

        return strtold(ch,(char **) NULL);
};

/*inline*/ double max(double a, double b)
{  return a>b ? a : b;
};

/*inline*/ double min(double a, double b)
{ return a<b ? a : b;
};

double dist(double A[3],double B[3])
{ long double res;
	
	res=sqrt(pow(A[0]-B[0],2)+
			 pow(A[1]-B[1],2)+
			 pow(A[2]-B[2],2));
	
	return res;
}
/////////////////////////////////////////histo
hist_type hist(float *ptr, int MAX_PTR, float lng[2]/*min, max*/,float delta){
	int i,j;
	hist_type sal;
		
	if(delta==0){
		if(lng[1]==0)/// por si les pasamos los max min manuelamente
		for(i=0,lng[1]=lng[0]=ptr[0];i<MAX_PTR;i++){
			if(ptr[i]<lng[0])lng[0]=ptr[i];			
			if(ptr[i]>lng[1])lng[1]=ptr[i];
		}	
		delta=(lng[1]-lng[0])/sqrtf((float)MAX_PTR);
	}
	sal.Ndat=(int) floor(sqrt((double)MAX_PTR));
	cout << "\nNumero de bines: "<< sal.Ndat << "\n";
//	getchar();
//	sal.dat = (float *) calloc((int)sal.Ndat,sizeof(float));
	sal.dat = new float[sal.Ndat+1];
	for(i=0;i<sal.Ndat;sal.dat[i++]=0);
	sal.max=lng[1];
	sal.min=lng[0];
	sal.delta=delta;
		
	for(i=0;i<MAX_PTR;i++){
		j=(int)floorf((float)sal.Ndat*(ptr[i]-sal.min)/(sal.max-sal.min));
		if(j>=0 && j<=sal.Ndat)
		sal.dat[j]++;
	}
//	cout << "salio!!!";
//	getchar();
	return sal;
}
int nside2npix(const int nside) {
  return (int) 12*nside*nside;
}
